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ABSTRACT 

The evolution of the metal content of galaxies and its relations to other global prop¬ 
erties [such as total stellar mass (M*), circular velocity, star formation rate (SFR), 
halo mass, etc.] provides important constraints on models of galaxy formation. Here 
we examine the evolution of metallicity scaling relations of simulated galaxies in the 
Galaxies-Intergalactic Medium Interaction Calculation suite of cosmological simula¬ 
tions. We make comparisons to observations of the correlation of gas-phase abundances 
with M* (the mass-metallicity relation, MZR), as well as with both M* and SFR or 
gas mass fraction (the so-called 3D fundamental metallicity relations, FMRs). The 
simulated galaxies follow the observed local MZR and FMRs over an order of mag¬ 
nitude in M*, but overpredict the metallicity of massive galaxies (log M* > 10.5), 
plausibly due to inefficient feedback in this regime. We discuss the origin of the MZR 
and FMRs in the context of galactic outflows and gas accretion. We examine the 
evolution of mass-metallicity relations defined using different elements that probe the 
three enrichment channels (SNII, SNIa, and AGB stars). Relations based on elements 
produced mainly by SNII evolve weakly, whereas those based on elements produced 
preferentially in SNIa/AGB exhibit stronger evolution, due to the longer timescales 
associated with these channels. Finally, we compare the relations of central and satel¬ 
lite galaxies, finding systematically higher metallicities for satellites, as observed. We 
show this is due to the removal of the metal poor gas reservoir that normally surrounds 
galaxies and acts to dilute their gas-phase metallicity (via cooling/accretion onto the 
disk), but is lost due to ram pressure stripping for satellites. 

Key words: cosmology: theory - galaxies: evolution - galaxies: abundances - galax¬ 
ies: haloes - galaxies: high-redsliift galaxies: star formation method: numerical 
simulations 


1 INTRODUCTION 


The study of elemental abundances in galaxies (‘metallic¬ 
ities’) can provide important information about their star 
formation histories as well as the different processes af¬ 
fecting gas evolution such as the ba l ance betwee n gas in- 
fall and outflows (e.g., Dave. Finlator fe Q ppenheimeijl201ll : 
IPaval, Ferrara fe Dunlodl2013l : lLillv et alj|2013h '. Moreover, 
the determination of correlations between the metal content 
of these systems and other properties such as their stellar 
masses (M*) or gas fractions (/ g ) is crucial for reconstruct- 
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ing the formation history of galaxy populations and for con¬ 
straining galaxy formation models. 


Galaxies in the local Universe, z ~ 0, exhibit a clear 
correlation between stellar mass and metallicity, such that 
systems with higher M, are more metal-enriched. This rela¬ 


tionship between A/* and metallicity (t he ‘mass-met alli city ’ 
relation, MZR) was first reported by iLeaueux et all (|l979l ) 
for star-forming galaxies, with subsequent observational 
studies typically us ing l u minosity as a surrogate (e.g., 


Garnett &; Shields 19871: IZaritskv, Kennicutt fe Huchral 
1994 Lamareille et al.ll2004 ). Over the last decade or so, im¬ 


proved data quality and larger statistical samples of galaxies 
have allowed for the deter mination of the MZR not only 
in the local Universe (e.g., iTremonti et al.ll2004 [Lee et al.l 
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20060 Kewlev & Ellison 20080) but also at higher redshifts 

up to 2 ~ 3.5 (e.g.. Savaglio et al. 

20051: Erb et al.l 20061: 

Maiolino et al. 20081: Havashi et al. 

20090: Lamareillc et al. 

20091: Mannucci et al. 2009: Lara-Lopez et al. 

201C; 

Moustakas et al. 2011: Cresci et al. 

20120: Yabe et al.l 

2012; 

Yuan, Kewlev & Richard 2013: Zahid et al.l 20141'). Q.ualita- 


tively speaking, these studies have typically found declining 
metallicities (at fixed stellar mass) with increasing redshift, 
although there are quantitative differences in the reported 
evolution rates. 

The scatter in the observed MZR has als o been 
extensiv ely studi ed rece ntly (e.g. jElhson_e£aJj 
Yates. Kauffmann fe Guol 120121 : IZahid et al 


alj 

200E 

. 

2012 


Hughes et alj 20131 ') as it encodes information about 


the dependence of metallicity on s e conda ry ph ysical param¬ 
eters. For example, ICooper et al.l ll2008ll and lEllison et all 
(1200911 have shown that, at fixed stellar mass, galaxies 
inhabiting denser environments tend to have higher metal¬ 
licities. Furthermore, higher gas-phase metallicities have 
been associated wi t h lower star formation rate s (SFRs, e.g., 
lEllison et al.l 20081: lAndrews fe Martinil 120130 ) and specific 


SFRs (e.g., 


Ellison et al 


2008), lowe r gas fractions (e.g. 


Zhang et al. 120091; Hugh es et all 201 30 ). smaller sizes (e.g., 
Ellison et all 20081: IWibe et al.l l2014l i . higher dust fractions 


(e.g., Zahid et al.l 20l3f and larger B — R and R — H 

air 


colours (e.g., Yabe et alll2O140 ). 

The correlation of the scatter in the MZR with other 
physical quantities suggests that the MZR is a 2D projection 
of a more fundamental 3 D (at least) relation. F or ex ample, 
lLara- Lopez et ahl (120100 ) and iMannucci et al.l (120100 ) have 
suggested the 3D relation betwee n M,, SFR and gas-p hase 
metallicity is more fundamental. [Mannucci et alj ( 20100 ) in 
particular analysed the metallicity as a function of M» and 
SFR. defining what they called the ‘Fundamental Metallicity 
Relation’ (FMR), in analogy wi th the fundamental p lane of 
elliptical galaxies introduced by [Dressier et all (1 19871 1. Since 
this 2D plane in 3D space is not parallel to the M* — Z plane, 
projecting the FMR onto the M* — Z plane introduces scat¬ 
ter. 


IMannucci et alj (l201oll defined a new quantity /x a = 
log(M*) — alog(SFR) with the aim of obtaining the pro¬ 
jection of the MZR with the minimum scatter, finding the 
tightest relation at a ~ 0.32. The FMR does not appear to 
evolve significantly up until 2 ~ 2.5, after which the metal¬ 
licity of galaxies of f ixed M* and S FR. sta rts to decline with 
increasing z (e.g., iTroncoso et al.l 120141 3 . Since the FMR. 
consists of a positive correlation between M* and metallic¬ 
ity at a given SFR and a negative correlation between SFR 
and metallicity at a given M,, part of the observed evolu¬ 
tion of the MZR. at 2 < 2.5 might be an artefact because 
galaxy populations that inhabit different regions of the FMR. 
are sampled at different cosmic epochs. For example, many 
high -2 surveys are magnitude-limited in the UV, tending to 
select high-SFR systems. As the metal content of galaxies 
decreases with SFR at a given A/*, at least part of the ob¬ 
served evolution of the MZR may be a consequence of the 
higher SFRs of galaxies obser ved at h igh 2 . Th is is consis¬ 
tent with the recent findings of lStott et al.l (l2013f j . who select 
systems with similar M* and SFR to local samples and find 
no significant evolution out of the MZR out to 2 ~ 1.5. 

While the existence of a 3D relation between M«, SFR 
and metallicity has been verified in recent years by many 


groups, who have explored different redshifts and mass 


of the relation (e.g 


Mannucci, Salvaterra & Camnisi 

2011; 

Cresci et al. 

2012 

Yabe et al. 20121: Cullen et al. 

2014; 

Henrv et al. 

2013! 

Lara-Lopez. Lopez-Sanchez & Ho 

pkins 

20ia 

20130; Lara-Lopez et al. 2013b 

; Pilyugin et al. 

Stott et al. 2013: Yabe et al. 2014 

). These discrepancies are 


generally related to different kinds of systematic uncer¬ 
tainties and selection biases that affect the observational 
studies and make the comparison between them non-trivial. 
For example, the use of different metallicity calibrations, 
methods for stellar mass determination and the aperture 
size may significantly influence the shape and zero point of 
the MZR. Moreover, at high 2 , the exploration of a wide 
parameter space is more challenging because of the limited 
sample sizes. The evolution of the physical conditions of 
star forming regions with cosmic time should also be taken 
into account as they can produce a differential evolution 
of metallicity indicators and locally calibrated metalli city 
relations may not be valid at high 2 llCuIlen et al.ll2O140 ). 

There is also still some debate about whether the M»— 
SF R-metallicity re l ation is truly fundamental. For exam ¬ 
ple. iBothwell et ahl (I2OI3I ') (see also lLara-Lopez et al.120133 ) 
proposed that the FMR can be considered a manifesta¬ 
tion of a more fundamental relation between A/*, HI mass 
and metallicity (the HI-FMR). These authors defined a new 
parametrization for projecting the MZR. in this 3D space: 
rip = log(AL*) — /31og(Af(HI) — 9.80), obtaining the lowest 
scatter at /3 = 0.35. While the dependence of the metallicity 
on SFR saturates at high masses, the anti-correlation with 
HI mass persists. Motivated by this work, we explore below 
3D trends between Ai», metallicity and both SFR and gas 
mass fractions. 

Various evolutionary scenarios have been proposed 
to explain the presence and evolution of the MZR (e.g. , 


Tisscra, De Rossi & Scannapieco 2005; Brooks et al. 2007; 


de Rossi. Tissera fe Scannapiecd 120071: iFinlator fe Davi 


2008 


DwAFinlafp r&^Oppenhameil 


2011 


Peeples fe Shank ar! 1201 ll: iDave. Finla tor fe Qppenheimeil 
2012|; iDaval. Ferrara fe Dunlopl [20Kf[ 


Yates. Kauffmann fe Guol 20121 : Romeo Velona et al.l 

2013ll . The lower metallicity of low-mass galaxies is 


generally thought to be a result of more efficient star 
formation-drive n outflows from shallo w er potentia l wells 
( Larsonj 197j; Tremonti et alj 20041 : iDalcantonl 120071 : 
iKobavashi, Springel fe White! 20070 ). However, a downsiz¬ 
ing scenario, with smaller galaxies exhibiting lower star 
formation efficiencies and consequently lower chemical en¬ 
richment, can also pot entially lead to a correlation between 
mass and metallici t y dBrooks et al.l 120071 : iMouhcine et al.l 
120081 : ICalura et al.l l2009h . Moreover, the infall of metal- 
poor gas from the IGM or through merger events has 
been frequently invoked as a key in gredient to repro¬ 
duce the observed trends ('e.g. iKoppen fe Edmundsl Il999l: 


Dalcanton. Yoachim fe Bernstein 120041: IFinlator fe Dava 

20081 : Dave. Finlator fe Qppenheimeil 201ll ~). In this sce- 
nario, the SF process is driven mainly by metal-poor 
inflowing gas, and the dilution of metal abundances by 
adding gas and the increased SFR. resulting from a bigger 
gas reservoir, then naturally leads to an anti-correlation 
between SFR and metal mass fraction Z. It is possible that 
differences in the morphological mix of galaxies and the 
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Metallicities of simulated galaxies 3 


associated changes in star formation ef ficiencies over time 
affect the evolution of the MZR (see e.g. ICalura et al.|[2009l . 
for a study in this sense). 


Aithough inflow and outflow of gas is likely to con¬ 
tribute to establishing the MZR, a complication to theoret¬ 
ical modelling is that the mass-loading of a wind, M w /M+, 
is not necessarily the same as the metal mass loading of a 
win d , Mz/M*. The su b-pc wind simulations presented by 
ICreasev. Theuns fe Bowen (l2015lf demonstrate that dense 
galactic disks drive metal enriched winds, whereas lower den¬ 
sity disks drive winds with far lower metallic ities. Interest¬ 
ingly, ee Rossi. Tissera fe Scannapiecd (120071 1 showed that 
a correlation between metallicity and mass naturally arises 
in a hierarchical framework as a consequence of the regu¬ 
lation of the SFR by merger events, without the necessity 
of invoking SN-driven outflows. However, without SN feed¬ 
back, they were not able to reproduce the observed slope of 
the relation (and neglect of SN feedback also leads to sig¬ 
nificant overcooling). Finally, the MZR need not result from 
the efficiency of winds as function of stellar mass, but might 
simply reflect variations in the stellar init ial mass function 
(IMF) (jKoppen, Weidner fe Kroupall2007l ). 


Although significant progress has been made both char¬ 
acterising and understanding the nature of the MZR and 
of the FMR, no clear consensus has emerged for either 
of them. In this paper we investigate if current models 
of galaxy formation and evolution are able to describe 
the observed trends, and if so, why. To do so we anal¬ 
yse the abundance evolution of star-forming regions in 
galaxies from the Gala xies-Intergalactic M edium Interac¬ 
tion Calculation (Mimic. [Crain et al.ll2009i ) suite of cosmo¬ 
logical hydrodynamical simulations. These simulations pro¬ 
vide a very large sample of galaxies allowing extrapola¬ 
tion of statistics to the (500 Mpc /i -1 ) 3 Millennium Sim¬ 
ulation dSpringel et al.l[20051 1 volume. They are able to re¬ 
produce rare systems such as voids and massive clusters 
and have been shown to describe well dynamical, kinemat- 
ical and abundance properties of obse r ved galaxies (e.g ., 
I Font et all 1201 ll : iMcCarthv et al.l 1 20121 : ISales et ail I2012T1 , 
and are therefore well suited to the problem at hand. We 
examine galaxy metallicity scaling relations and their evolu¬ 
tion in the GIMIC simulations. In particular, we make com¬ 
parisons between GIMIC and observations of the MZR (and 
its 3D extensions) of local galaxies, we explore the evolution 
of metallicity scaling relations with stellar mass and circular 
velocity, and we examine the difference between central and 
satellite galaxies in terms of their enrichment histories. We 
also analyse the abundance evolution of different elements. 


The plan of the paper is as follows. The simulation and 
sample of galaxies are described in Sections 12.11 and 12.21 
respectively. The redshift z = 0 mass-metallicity relation is 
presented in Section l3TTl and compared to observations, while 
in Section 2] we analyse and discuss its evolution. In Section 
[5] we compare the metal enrichment of central and satellite 
galaxies in GIMIC. In Section [Gl we examine the dominant 
processes that determine the main features of the simulated 
mass-metallicity relation. Finally, our conclusions are sum¬ 
marised in Section [7] 


2 THE gimic SIMULATIONS 


2.1 Cosmological parameters and sub-grid 
implementation 

ICrain et al.l d2009i ) provides full details on the GIMIC sim¬ 
ulations, here we begin by giving a brief overview. Five 
nearly spherical regions of radius R ~ 20 Mpc /i _1 where 
picked at redshift z = 1.5 from the Millennium simulation 
JSpringel et al.l[2005l l. and resimulated using ‘zoomed’ initial 
conditions and including hydrodynamics. In such zoomed 
initial conditions, the picked regions are simulated at high 
resolution, with the rest of the (500 Mpc h -1 ) 3 Millennium 
volume represented much more coarsely with more massive 
dark matter particles. The five regions were selected to have 
mean overdensities that deviate by (-2, -1, 0, +1, +2) a 
from the cosmic mean, where a is the root mean square 
mass fluctuation on a scale of ~ 20 Mpc /i _1 at z = 1.5. 
These five regions sample environmentally diverse cosmo¬ 
logical regions, ranging from massive clusters to voids. Such 
a zoom technique allows us to combine the advantages of 
having high numerical resolution with sampling a range of 
cosmological diverse environments. 

The simulations assumes a cold dark matter uni¬ 
verse at present dominated by a cosmological constant, 
and uses the same values for the cosmological parame¬ 
ters as the Millennium Simulation, (fl m , Ha, ug, n s ) = 
(0.25, 0.75, 0.045, 0.9,1), where the symbols have their usual 
meaning. The simulations were performed with the code 
ga dget-3, a mod ified version of gadget-2 last described 
by 1 pringell (|2005l f. and with the sub-grid ingredients for 
galaxy formation taken fro m the O verwhelm ingly Large 
Simulations (OWLS) suite dSchave et al.l l2Qloh as summa¬ 
rized below. Hydrodynamics in GADGET code is based on the 
smoot hed particle hydrodynamics (SPH ) scheme of iLucvl 
d 197711 and iGingoId fe Monaghanl (1 19771) . with the cosmic 
gas represented by a set of particles. 

Radiative cooling rates for cosmic gas are com- 
puted element-by-elem ent in the presence of a 
iHaardt fe Madaij (120011 1 ionising UV/X-Ray background 
as a function of density, temperatur e and redshift us¬ 
ing C LOU DY iFerland e t _alj ll998l f as described in 
IWiersma, Schave fe Smithl ( 20091 1. Reionisation of Hydro¬ 
gen is assumed to take place at redshift z = 9, and of 
Hell at redshift z = 3.5, motivated b y observ ations of the 
intergalactic medium (e.g.. iTheuns et alll2002r i. 

St ar formation is implemented fo llowing the prescrip¬ 
tion of ISchave fc Dalla Vecchial (120081 ). Above a hydrogen 
number density threshold of p > 0 .1cm~ 3 , gas can self-shield 


from the UV-background (e^g., Altav et al . 201 il l . and be¬ 
comes thermally unstable dSchavel 2004l f. Lack of numeri¬ 
cal resolution prevents the GIMIC simulations from follow¬ 
ing this process in detail, and our sub-grid implementation 
parametrises the unresolved physics by increasing the pres¬ 
sure of such star forming gas to p = po (p/po) 4//3 . Taking the 
star formation rate to be the power-law 


P* = P*,o (p/p 0 )" , 


(1) 


yields simula ted g alaxies th at follow the Kennicutt- 
Schmidt law (IKennicut tl 119981 ) when in hydrostatic equi¬ 
librium, with the constants p+,o and n directly related 
to the amplitude and slope o f the Kennicutt-Schmidt law 
(ISchave fe Dalla Vecchial 12008 1. Star formation is then im- 
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plemented stochastically by converting SPH particles to col¬ 
lisionless ‘star’ particles with a probability set by their star 
formation rate. Note that, in the case of the very recent EA¬ 
GLE simulations, SFR is implemented in a very similar way 
b ut cons ideri ng a metallicity dependent SF thresdhold (see 
ISchave et alj|2015l . for details). 

Abunda nce evolution an d enric hment is implemented as 
described by IWiersma et all d2009h . The simulation tracks 
eleven elements (H, He, C, Ca, N, O, Ne, Mg, S, Si, Fe) 
Q, produced by three stellar evolutionary channels, (i) core¬ 
collapse (type II) supernovae resulting from massive stars 
(M > 6Mg), (ii) type I supernovae assumed to result from 
catastrophic mass transfer in close binary stars, and (in) 
asymptotic giant branch (AGB) stars. Stellar evolutionary 
tracks and yields_that depend on initial me tal ab undance 
are ta ken fr om IPo rtinari. Chiosi fe Bres san ( 19981): iMarigol 
ll200lll ; and iThielemann et al.l 1 20031 ). 

Nucleosyn thetic yields are uncert ain by factors of a few 
as discussed in lWiersma et aD J20091 ) (see also the detailed 
comparison of yield s for a range o f elements in published 
models by iNomoto, Kobavashi fe Tominagal 120131 ) . Apart 
from uncertainties in the yields of a given element for a 
given stellar evolutionary channel, there is also discussion in 
the literature about which channels domina te for a given el¬ 
ement. For example. Ichiappini et all J 2 OO 6 I ) argue that rota¬ 
tion at low metallicity significantly enhances the production 
of Nitrogen and that this is required to solve the primary Ni¬ 
trogen prob l em for low-metallicity Milky Way stars whereas 
iKobavashi fe Nakasatd (l201lh argue that AGB stars could 
b£_the source of t he N even at low met allicity as discussed in 
iNomoto, Kobavashi fe Tominagal (j2013l ). A detailed investi¬ 
gation of the abundances of different elements in different 
en vironment s could help to resolve such issues. For example 
IPipino et al.l (120091) shows how the changes in N yields as a 
function of metallicity that resolve the primary N problem 
in the Milky Way fails to predict the correct N/Fe ratio in 
local spheroids. In this paper we use a single set of yield 
tables and show that abundance patterns of galaxies never¬ 
theless vary significantly between them because the extent 
to which a galaxy can hang-on to its own metals depends 
on the efficiency of feedback. This clearly complicates mat¬ 
ters and illustrates how one cannot simply compare yield 
tables to observations and judge which table is the more ac¬ 
curate. Changing yield tables in addition to making changes 
in the efficiency of feedback will of course also impact the 
abundance evolution. 

A star particle in the simulation represents a single stel¬ 
lar p opulatio n with given age and initial abundance, assum¬ 
ing a lChabrieij (120031) IMF. Star particles distribute the syn¬ 
thesised elements and total mass lost during each time-step 
to neighbouring gas particles according to the SPH interpo¬ 
lation scheme. 

Kin etic f e edback from SNe is implemented as described 
bv lDalla Vecchia fe Schavel (120081 ). After a time delay of 3 x 
10 7 yr, the mean life time of a core collapse progenitor, a 
star particle ‘kicks’ on average 17 = 4 gas particles into a 
wind with launch speed v w = 600 km s _1 . The values yield a 
star formation history that tracks the observations well. The 


1 [Wiersma, Schave fc Smithl (2009) found that these elements 
contribute significantly to the radiative cooling at T > 10 4 K. 


chosen values for the parameters ( 17 , v w ) imply that 80 per 
cent o f the avai lable SN energy is used to drive a wind (for a 
IChabrierl (120031 ) IMF with stars forming with masses in the 
range 0.1-100 Mg, of which those more massive than 6 Mg 
undergo core collapse). Energy injected by type I SNe is also 
taken into account. It is worth noting that in this scheme, 
the ’kicked’ particles are not hydrodynamically decoupled 
at any stage. 

The GIMIC simulations have been run at three levels 
of resolution. The low-resolutio n versi o n correspo nds to the 
original Millennium simulation (ISpringel et al.l2005l) : the in¬ 
termediate and high resolution correspond to gas particle 
masses of 1.16 x 10 7 and 1.45 x 10 6 /i~ 1 Mq, respectively, 
with the dark matter particles a factor of — fib/flb = 4.56 
higher. Gravitational forces are spline-softened with equiv¬ 
alent Plummer-sphere softening length of t = 0.5/i _1 co¬ 
moving kiloparsecs for the high resolution run until redshift 
z = 3, and is kept constant in physical units thereafter. The 
intermediate resolution simulation uses e = l/i _1 kpc. 

For most of the analysis presented below we used the 
high resolution GIMIC simulation performed in the -2 and 
Ocr regions. The Appendix compares these results with those 
from the intermediate resolution runs; we find that the main 
trends and conclusions are robust against a change in nu¬ 
merical resolution of a factor of eight in mass. Below and 
unless otherwise stated, we will refer to the mass fraction of 
elements heavier than Helium as ‘metallicity’. 


2.2 Identifying galaxies in GIMIC 

We identify dark matter halos by applying a friends- 
of-friends (FoF) algorithm with linking length b = 
0.2 times the mean interparticle separation dDavis et alJ 
1985 ). We then apply the s ubfin d algorithm described in 


Springel . Yoshida fe W hitel (l 200 ll) and modified to include 
barvons by iDolag et al.l 1 20091) to identify galaxies as self¬ 


bound structures. We will refer to the most massive SUB- 
find structure of a FoF halo as its ‘central galaxy’ , any 
other galaxy in the same FoF halo as a ‘satellite’. Unless 
otherwise specified, central substructures of FoF-halos are 
analysed in this work. To make comparison of simulated 
galaxies to observations more realistic, we measure prop¬ 
erties of simulated galaxies out to a maximum distance of 
Rcut = 20 kpc to mimic surface brightness limitations in the 
data. 

We limit our analysis to subst ructures wit h mo re than 
2000 particles in total, and following lMcCarthv et al~](l2012l) . 
also require galaxies to have a stellar mass of at least 
M* 10 9 AIq (approximately 500 star particles). We also 
require M* ^ 10 1 o ' 5 Mq, because GIMIC does not include 
AGN feedback and as a consequence these more massive 
galaxies are overcooled and have too high stellar fractions. 
These selection cuts result in a final sample of between 450 
and 650 simulated galaxies, depending on redshift. We ver¬ 
ified by using the 5 spheres of the intermediate resolution 
GIMIC suite that our findings do not depend (strongly) on 
large-scale-environment and hereafter we do not distinguish 
between the two different high-resolution GIMIC spheres (i.e. 
the -2 and Oct GIMIC runs). 

Finally, element abundances are stored based on both 
particle and SPH s moothed quantities as explained in 
IWiersma et al.l (|2009l ). The SPH algorithm computes the 
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density p(r^) at the location ri of particle i by summing 
the kernel W over all particles within a smoothing length hi 
of 

p(u) = W(\rirj\/hi ), (2) 

3 

where rrij is the mass of particle j. If Aj is the mass fraction 
in a given element for particle j (i. e. its ‘particle metal- 
licity’), then the smoothed metallicity Ai for element A of 
particle i is calculated as 

AiPin) = 'Y^A j rn j W{\TiTj\/hi). (3) 

3 

Our results do not depend on which of the two abun¬ 
dances we use, for simplicity we only show results corre¬ 
sponding to smoothed abundances below. 


3 COMPARISON TO OBSERVED 

ABUNDANCE RELATIONS 

3.1 The mass-metallicity relation 

The Af*-Z mass metallicity relation (MZR) is a fundamen¬ 
tal galaxy scaling relation, linking dynamics in terms of star 
formation, gas inflow and outflow, to abundances. In this 
section we compare the MZR in GIMIC to the observed rela¬ 
tion at redshift z = 0, which is typically inferred by measur¬ 
ing abundances in star forming gas. We will therefore begin 
by looking at gas-phase abundance in GIMIC weighing each 
particle with its star formation rate computed as described 
in Section [2] The SFR weighted abundance A of an element 
A of a galaxy is 

A Af* = ^ , (4) 

i 

where the sum is over all gas particles i in that galaxy, rh*,i 
and Ai are the star formation rate and abundance of ele¬ 
ment A of gas particle i, and AT* = ^. m*^ is the total star 
formation rate. For simulated galaxies we calculate metallic¬ 
ities within the inner 20kpc, but verified that using 10 kpc 
instead makes no significant difference. 

GIMIC galaxies exhibit a tight correlation between gas- 
phase Oxygen abundance weighted by star formation rate, 
and stellar mass, with more massive systems more metal rich 
(Fig. [T]). The ratio O/H increases by approximately 1 dex 
over the mass range 10 9 — 10 10 ' 5 Mq, with approximately 

1 dex of scatter at given A/*. _ 

_ The observed values of O/H quoted by[Mannucci_et_£A 

(2010|) are ~ 0.4 dex higher than those of ISanchez et all 
(2013), but the dependence on Af* is very similar, with O/H 
increasing with Af* up to Af* ~ 10 10 ' 5 Mq, becoming ap¬ 
proximately constant_atJii£her mass. The O/H-A/* relation 
reported bv lLara-Lopez et al.l (l2010l l is sign ificantly steeper. 
GIMIC metallicities follow the O/H-Af* from lMannucci et al.l 
(l2010ll up to Af* ~ 10 10 AfQ, but are higher than those 
data for more massive galaxies. The absence of AGN feed- 
back in GIMIC causes overcooling in more massive objects 
l|McCarthv et al.ll2012f) . and the resulting high stellar frac¬ 
tions and lack of feedback may be responsible for the high 
O/H values. 


AGN-feedback plays a crucial role in quenching star 
formation in massive galaxies (e.g. iBower et al]|2006ll and 
models that do not consider AGN tend to over predict st ellar 
mass and star formation rates in such galaxies. ISchave et al.l 
d2015l l present results from the EAGLE0 simulations which 
reproduce the observed galaxy stellar mass function well and 
show that these simulations agree well with the observed 
MZR also for more massive galaxies (their Fig 13). We plan 
to investigate the metal-enrichment history in EAGLE in a 
future work. 


3.2 O/H dependence on SFR and gas fraction 

At a given stellar mass, the O/H abundance of GIMIC galax¬ 
ies decreases with increasing gas fraction, f g = M g /(M * + 
M g ) = Mg/Mb, and increasing star formation rate (Fig. [l] 
middle and right hand panels, respectively). Here, f g refers 
to the total gas fraction - not just the fraction of star forming 
gas - but we verified that similar trends appear when lim¬ 
iting the gas mass to star forming gas only for lower mass 
galaxies with Af* < 10 10 ' 3 Afg. The limited sample size re¬ 
stricts us to using a small number of bins in f g and SFR, 
and each bin contains 30-40 per cent of the full sample. At 
the low-mass end of the simulated MZR, the difference in 
O/H between the highest and lowest f s and SFR bins are 
~ 0.3 dex and ~ 0.2 dex, respectively. 

The decrease in O/H with increasing f„ and SFR 


is also seen in the data of iMannucci et al. (l20ld l and 
iBothwell et al.l (l2013h . The solid lines in the right h and panel 
correspond to the FMR of IMannucci et all (l2010lh with red 
and pink lines corresponding to log(SFR) = —2 and 0, re¬ 
spectively. Encouragingly, we find reasonable agreement in 
zero-point and slope of the simulated relations compared to 
data up to Af* ~ 10 w Mq. At higher Af* ^ 10 10 ’ 5 Afg GIMIC 
galaxies are too metal rich as noted previously. 

In GIMIC, SFR is related to the gas reservoir through 
the imposed (subgrid) star formation law of Eq. |l] see 
Section E 3 >, which does not simply translate to Af* oc M g . 
However M g and SFR are correlated in GIMIC galaxies, with 
galaxies with a larger total gas mass exhibiting higher SFR. 
By itself, this does not imply a causal connection: the cor¬ 
relation may be a consequence of the same phenomenon, for 
example an increased gas fraction and enhanced star forma¬ 
tion may both result from a recent accretion or merger event 
that adds gas to the halo. Whatever the origin of the gas, 
some of it has low elemental abundance and makes it into 
the galaxy, since the abundance of star forming gas decreases 
with increasing SFR. This suggests that at lower stellar mass 
where these correlations are strongest, accreted low abun¬ 
dance gas causes both the decrease in abundance, and the 
increase in star formation. The good agreement between the 
scaling relation in GIMIC and the observations then suggests 
that the same process may be at work in real galaxies. 

Interestingly, the observed relations turn-over for more 


2 EAGLE (Evolution and Assembly of GaLaxies and their Envi¬ 
ronments) is a project of the Virgo consortium for cosmological 
supercomputer simulations. More information can be found at the 
EAGLE web sites at Leiden, http://eagle.strw.leidenuniv.nl/, and 
Durham, http://icc.dur.ac.uk/Eagle/ 
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Figure 1. SFR-weighted gas-phase Oxygen abundance, O/H, as a function of galaxy stellar mass, M*, at redshift z = 0. Left panel: 
full GIMIC sample of central galaxies (small symbols) with mean and 25 and 75t h percentiles shown as b l ack solid line with error bars. 
The red solid , blue d ashed and pink dot-dashed lines are observational fits from [Mannucci et alj (120101) ; lLara-Lopez et ab I tepid) and 
1 Sanchez et al.l (I2013T) , respectively. There are significant systematic differences between the observational fits, both in amplitude and 
shape. The simulated abundances are consistent with the data at lower masses but may be slightly higher at the most massive end, 
M* « 10 10-5 Mq. Middle panel: Same as left panel but for GIMIC galaxies binned in gas fraction, f g = M g /(M g + M*) = M g /M At 
given M*, gas rich galaxies have lower O/H. Right panel: same as left panel, but for ga laxies binned in star fo rmation rate, M*. At given 
M*, galaxies with higher M* have lower O/H, consistent with the observed trend from [Mannucci et al.l (120101) shown for log(SFR) = —2 
and 0 as red and pink line, respectively. 



Figure 2. Location of observed galaxies in M* - SFR - O/H space from [Mannucci et alj (120101) ( left panel) and of GIMIC galaxies (righ t 
panel). The shaded area in the left panel represents the best fit 2D surface to a fundamental plane taken from [Mannucci et al.1 (j2010h . 
the right panel gives our fit to GIMIC galaxies. Lines with symbols depict the median M*-0/H relations for bins in SFR from Fig. [T] (right 
panel), black lines trace the shaded surface at the center of each SFR bin. The relation between M* and SFR in the simulations is shown 
by the symbols in the M*-SFR plane. 
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massive galaxies, M ~ 10 10,5 AfQ, in the sense that metallic- 
ity no longer increases with A/*. As gas surely still accretes 
onto the haloes of these more massive galaxies, the break 
down of the scaling relation then suggest that the gas no 
longer accretes onto the galaxy. Contrasting these observed 
relations to what happens in GIMIC where the scaling rela¬ 
tion does hold for these more massive galaxies suggests the 
cause of this: feedback. In real galaxies some process pre¬ 
vents the gas from accreting onto the galaxies, but not in 
GIMIC. We suggest it is likely AGN feedback - not included 
in GIMIC - that causes the turn-over in the observed MZR. 

The location of observed galaxies in O/H, A/* and 
SFR. space is plotted in Fig. [2] ( left panel). The correla¬ 
tions discussed so far cause observed galaxies to scatter 
ar ound a (curved) 2D plane in this 3D space, and the fit 
of lMannucci et ahl (l2010li to this ‘fundamental plane’ 

12 + log (O/H) = 8.90 + 0.37m - 0.14s - 0.19m 2 

+ 0.12ms- 0.054s 2 , (5) 

is shown in the left panel. Here, m = log(M*/10 10 Mq) and 
s = log(SFR/MQ yr _1 ). Similar to the observed galaxies, 
GIMIC galaxies also trace a plane in this space, with 

12 + log (O/H) = 8.93 + 0.79m - 0.047s + 0.0048m 2 

+ 0.19ms-0.015s 2 . (6) 

The ratio O/H depends more strongly on mass and less 
strongly on star formation rate in the simulations as com¬ 
pared to the data, but the overall similarity is encouraging. 

This GIMIC fundamental metallicity relation describes 
well the median A/*-O/H relation at different SFRs. For 
both observed and GIMIC FMRs, O/H is anti-correlated with 
SFR at the low-mass end (M* < 10 9 Mq), but this trend re¬ 
verses for more massive galaxies (M* > 1O 1O M0) for which 
O/H is correlate d with SFR. Such trends ar e consi stent with 
those found by iYates. Kauffmann fc Guol (120121 ) in semi- 
analytic model galaxies. 

In a nalogy with the fu ndamental plane of elliptical 
galaxies (IDressIer et al]|l987i f. it is possible to chose combi¬ 
nations of O/H, SFR. and Af* that minimise the scatter when 
seeing the FMR edge on, for example using the combination 
O/H and p. a = log (A/*) — a log(SFR). For GIMIC, scatter is 
minimised for a v alue of a = 0.2 as com pared to a = 0.32 
for the data from iMannucci et al.l (l20ld ). Alternatively in 
terms of rjp = log(M*) — /31og(M(HI) — 9.80), GIMIC galax¬ 
ies prefer /3 tts 0.2 when using the mass in star forming gas in 
lieu of Af(HI), and /? = 0.48 when using the total gas mass, 
as compared to 0 = 0. 35 for observed galaxies according 
to lBothwell et all ll2013h . We conclude that the feedback in 
GIMIC induces similar correlations between abundance and 
star formation rate as a function of galaxy stellar mass, al¬ 
though in detail there are differences, especially for more 
massive galaxies. 


3.3 The / g — Af* - SFR relation 

The anti-correlation of the O/H ratio with gas fraction 
(/ g ) and SFR. at given galaxy stellar mass can be exam¬ 
ined in more detail by analysing the / g — Af* — SFR rela- 
tionship in GIM I C. To compare with observations we follow 
iDe Rossi et al.l (|2013l ) and select ‘warm’ gas with temper- 



9.0 9.5 10.0 10.5 

log(M.) [M 0 ] 


Figure 3. Median warm gas fraction, / g = A7 g (T < 
15000 K) /My, , as a function of M* at redshift 2 = 0, A/), is 
the galaxy’s baryon mass ( black solid line), with error bars in¬ 
dicating the 25th and 75th percentiles. The trend is in good 
agreement with the observed fraction from lStewart et all (2009) 
(dashed black line). Coloured lines show the median of / g in the 
simulation for bins in star formation rate, see legend. At given 
stellar mass, galaxies with higher gas fractions have higher SFR. 

aturt@ T ^ 15000 K, but note that the trends we find re¬ 
main almost the same if we use the total gas component. 
For ~ 85% of our sample, this ‘warm’ gas component rep¬ 
resents more than ~ 75% of the total gas content. Gas with 
T > 15000 K dominates in mass in only ~ 7% of GIMIC 
galaxies. 

The fraction of warm gas in the simulations falls with 
i ncreas ing stella r ma ss, and follows the observed relation of 
IStewart et al.l J2009I) very closely (Fig. [3j. At given M*, the 
simulated / g increases with increasing star formation rate, 
and galaxies with little or no sta r formation have very lit¬ 
tle warm gas. ISantini et al.l (120141 ) reported the existence of 
a fundamental / g — A/* — SFR relation in observed galax¬ 
ies, and the GIMIC galaxies appear to be consistent with 
this. Nevertheless, in the SFR bin where observations and 
simulations are both available, GIMIC galaxies tend to over 
predict gas fractions for massive galaxies. We discuss stellar 
metallicities next. 

3.4 Stellar metallicities 

The dependence of stellar metalli city on ste llar mass in 
GIMIC follows the data from iGallazzi et al.l (120051 ) very 
closely (Fig. [4] left panel), except for the most massive galax¬ 
ies (A/* > 1O 1O,S M0) which are overcooled in the simulation 

3 We reiterate that GIMIC imposes a pressure-density relation for 
dense gas, and hence does not distinguish between the warm and 
cold neutral medium in the ISM of a galaxy. 
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Figure 4. Median stellar metallicity as function of stellar mass (left panel), and as function of SFR weighted gas O/H (right panel). 
GIMIC metallicities are shown as symbols, with black line showing the median relation, and the 25th and 75th percentiles indicated as 
error bars. Observational data from Gallazzi et al. (2005) are shown as full red lines (median), with dashed lines encompassing the 25th 
and 75th percentiles. The simulated relation follows the observed relation well in the left panel, but the most massive simulated galaxies 
are more metal rich as noted before; the scatter in the data is larger than in the simulation. The observed relation in the right panel 
is steeper than in the simulation, with simulated galaxies reaching a plateau of log (Z+/Z@) ~ —0.4 at low gas abundance, which is not 
apparent in the data. 


0. Also, the scatter around the median is significantly larger 
in the observations. 

We have verified by measuring GIMIC stellar abundances 
in apertures [i.e. varying the value of the parameter R cu t 
discussed earlier), that the smaller scatter in GIMIC is not 
likely a consequence of aperture limits. However uncertain¬ 
ties in observed measurements of Z k may cause the intrinsic 
scatter to be overestimated (see IGallazzi et al.l 120051 ). In¬ 
cluding satellites in our sample does increase the scatter in 
the M* — Zi, relation, and we also find that aperture effects 
do affect satellites’ metallicities. It may be that observa¬ 
tional uncertainties, and the inclusion of satellites, would be 
enough to reconcile observe and simulated scatter, but we 
have not verified this in detail. 

In Fig. U ( right panel) we explore the relation between 
stellar metallicity Z+ and gas-phase O/H weighted by star 
fo rmation rate. We find reasonable agreement with the data 
of IGallazzi et al.l (120051) . with the possible exception of the 
behaviour at low abundance: the simulated Z+ reaches a 
plateau in Z* at low gas-phase O/H not apparent in the 
data. This trend in the simulations is consistent with a SF 
triggered by infall of metal poor gas: in Fig. [Tjwe saw that 
most metal-poor systems have the highest gas fractions, f s , 
and correspondingly higher star formation rates. 


4 METALLICITY EVOLUTION 

Previously we demonstrated that GIMIC galaxies display sim¬ 
ilar trends in ISM and stellar metallicity at redshift 2 = 0 
as observed galaxies. This reasonable agreement encourages 

4 Following IWiersma et all (2009 ). for the solar abundance, we 
use the metal mass fraction Zq = 0.0127. 


us to explore the abundance enrichment history. In previ¬ 
ous sections we have analysed only O/H as Oxygen is the 
metal most commonly used in observational studies. In the 
following, we will take advantage of the enrichment evolu¬ 
tion of all elements tracked in the simulation. We caution 
the reader that our results of course apply to the particular 
set of yields discussed in Section [2] 


4.1 Abundance evolution for different elements 

The abundance evolution of different elements in the simu¬ 
lation is plotted in Fig. [5] depicting gas phase abundances 
weighted by star formation rate. The O/H-M* relation ( up¬ 
per left panel) - commonly used as a proxy for the ‘mass- 
metallicity relation’ does not evolve in GIMIC either in slope 
or normalisation, and is in place at least since 2 ~ 3. We 
verified that this result does not depend on the aperture 
within which O/H is measured by changing the value of 
Rcut defined earlier. In addition, we employed the redshifts 
used in the observational an alysis of the evolution of the 
MZR of lMaiolino et all d2008h . estimating O/H inside their 
observed radii: 3.0 kpc (for 2 « 3.5), 3.6 kpc (z ~ 2.25), 
3.725 kpc (2 « 0.74) and 2.0 kpc (2 ss 0.06) to enable a 
fair comparison. Using these apertures has no effect on the 
slope of the O/H-M* relation, but the zero point increases 
by 0.2-0.3 dex - a consequence of metallicity gradients in the 
simulated galaxies. Including satellites in the analysis does 
introduce modest evolution in O/H with 2 , with higher 2 
galaxies having lower O/H at given stellar mass, see also 
Section [5] below. 

The absence of evolution in the MZR in GIMIC appears 
inconsistent with the results of several observational studies 
which claim significant evolution since 2 = 3, in the sense 
of higher redshift galaxies having lower O/H at given stel- 
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Figure 5. Median SFR-weighted gas abundances as a function of M* for different redshifts z (different line styles) for different elements 
(different panels); error bars depict the 25th and 75th percentiles of galaxies. To ease comparison between different elements, abundances 
are off-set by a constant - Y - whose value appears in brackets in each panel (Y = 12 for Oxygen - the usual normalisation for that 
element). Top panels are for a-elements produced in the simulations mainly by type II SNe, while bottom panels are for N and C, 
produced by AGB stars, and Fe produced in both type I and type II SNe. 


lar mass. For example. iMaiolino et all (1 2008 ) measure a de¬ 
crease in O/H for galaxies with M* = 10 10 Mq of 0.8 dex 
between redshifts z = 0.7 and z = 3.5 (their Figs. 7 and 8), 
with lower M* galaxies evolving even faster; at that stellar 
mass, z = 2.2 galaxies are less enriched than z = 0 galaxies 
by 0.4 dex. 

A plausible explanation for the lack of strong evo- 
lution in the MZR in the simulations can be found in 
IWeinmann et all (120121 ). These authors showed that many 
current cosmological simulations (including GiMic), as well 
as semi-analytic models, form too much stellar mass too 
early compared to observations. As the simulations have 
approximately the correct amount of stellar mass (per unit 
halo mass) at 2 ~ 0.1, this implies that the simulations will 
also underpredict the amount of late-time star formation. 
A likely consequence of this is that the MZR will be set at 
higher redshift in the simulations compared to real galaxies. 

Alternatively, or perhaps in conjunction with the above, 
the presence of a strong correlation between O/H-M* and 
M* in simulated - and observed - 2 ~ 0 galaxies suggests a 
possible bias plausibly present in the observations: if higher 
2 galaxies had higher M* due to observational selection, then 
their lower metallicities might be inferred to imply evolution. 


Indeed, observations by lStott et al.i (|2013i ) suggest a negli¬ 
gible evolution of the MZR from 2 = 0to2~l - 1.5, and 
these authors explain the discrepancy with previous studies 
due to selection bias. Such an explanation is consistent with 
GIMIC in which low-mass galaxies follow the FMR at least 
up to 2 = 2. 

Figure [5] also displays the evolution of the abundances 
of other nuclei (Xi). As expected Oxygen is the dominant 
metal by mass at all 2 and as such is a good proxy for 
metallicity. The shape of the Xi/H-M* relation and its evo¬ 
lution (or rather lack thereof) is almost identical for Oxygen 
and Magnesium, both of which are a-elements synthesised 
in massive stars and released by Type II SNe. However, 
the Si/H relation appears to flatten toward 2 = 0, even 
though Si is yet another a-element. In the yield tables that 
GIMIC uses for Type II enric hment (basically the yields from 
IPortinari. Chiosi fe Bressanl 119981 with some modifications 
suggested by L. Por tinari in a private communication, see 
IWiersma et al]|2009l for more details), the metallicity depen¬ 
dence of yields is similar for O and Mg: at higher metallicity, 
massive stars (M > 30 Mq) produce less O and Mg than at 
lower metallicity - whereas the behaviour for Si is the op¬ 
posite (i.e. at higher metallicity stars tend produce more 
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Figure 6. Median SFR-weighted gas-phase metallicity as a func¬ 
tion of M* for different redshifts; error bars depict the 25th and 
75th percentiles. At a given stellar mass, the ISM is higher en¬ 
riched at lower z although the trend is weak given the scatter, 
and is nearly absent for M* > 10 10 Mq . 


Si than at lower metallicity). However in all, there is little 
if any evolution in the mass-metallicity relation for any of 
the a elements O, Mg, and Si, with scatter in abundance 
significantly larger than net evolutionary trends. 

The bottom left panels of Fig.[5]show evolution of Nitro¬ 
gen and Carbon. In our models, AGB stars contribute signif¬ 
icantly to C, and they are the dominant source of N. Abun¬ 
dances of N/H and C/H show significant evolution particu¬ 
larly from z = 3 —> 1, with an increase in abundance towards 
lower z, accompanied by a flattening of the C/H-M* rela¬ 
tion . In the nucleosynthesis models we employ, the N yield 
is more sensitive to metallicity than the C yield at higher 
abundance. This is due to the secondary element nature of 
N - synthesised from other metals rather just from burning 
Hydrogen. This might cause the N abundance turn-up for 
more massive - and hence higher Z galaxies, with approx¬ 
imately N/H proportional to Z 2 . Recent chan ges to Nitro¬ 
gen y ields at low metallicity discussed by e.g. IPipino et al] 
1: 20091 ) could plausibly change the evolution of N and a more 
detailed comparison of the present models to data could test 
such a model. This is however, beyond the scope of this pa¬ 
per. 

The strongest evolution is for Iron (bottom right panel of 
Fig-0. and happens mostly below z = 2. Iron has a signifi¬ 
cant contribution from Type I SNe which release synthesised 
Iron with significant delay of ~ 1 Gyr after the formation 
of the progenitor star. At M* = 10 10 Mq, the Fe abundance 
increases by almost 0.3 dex between z = 3 and z = 0, and 
the increase is even more pronounced at lower M*. 

Of course these abundances are not just related to nu- 
cleosynthetic yields and how they depend on Z, but also to 
the ages of the population, and to what fraction of the el¬ 


ements synthesised in stars manage to remain in the ISM. 
GIMIC galaxies have on average similar mass-weighted stellar 
ages of ~ 10 Gyr, therefore comparing galaxies of a given 
stellar mass, those at lower 3 had more time to become en¬ 
riched with elements synthesised in SNIa and AGB stars, 
driving the observed evolution of non-o elements. The (unre¬ 
alistically) high ages of simulated galaxies are a consequence 
of too early star formation in these systems, and leave them 
with a significant old stellar component. The absence of a 
M*-age relation in the simulation might artificially enhance 
the correlations between abundances and redshift, and may 
also lead to an underestimate of the scatter in simulated 
abundances. 

The metallicity evolution of the gas phase, weighted by 
star formation rate, shows only very mild evolution, with Z 
increasing with time (Fig.0, with a scatter between galaxies 
of given Af* at any given time as large as the evolution of the 
median Z from 2 = 3 to 3 = 0. Similarly to the case of stellar 
abundances this is a consequence of the absence of signifi¬ 
cant evolution of O/H, and other elements may well show 
evolution. The different behaviour of Oxygen compared to 
several other elements suggest that caution should be exer¬ 
cised when using this element just by itself as a measure of 
abundance evolution. 


4.2 Metallicity as a function of circular velocity 

One theory of the origin of the MZR is that systems 
with lower Af* are m ore affected by galactic winds (e.g. 
iFinlator &; Davj l2()()8l ). Such ‘leaky’ galaxies cannot effi¬ 
ciently retain their ISM gas or recently produced met¬ 
als, with SNe ejecting both gas and metals from the disk 
or indeed the halo of such galaxies (see discussion in 
ICreasev. Theuns fc Bowen 120151 '). Such efficient feedback is 
a crucial ingredie nt in models of galaxy evolution (e.g., 
ISchave et al.|[20lol l , and GIMIC includes a sub-grid model for 
launching winds. The extent to which gas mass is ejected 
from a galaxy by winds is thus more likely to depend on the 
circular velocity of the galaxy - a measure of the depth of 
the potential well - rather than its stellar mass, although 
the gas density and hence the fraction of energy in the 
wind lost through radiative cooli ng may play a role as well 
dCreasev, Theuns fe Bowerll2015h . 

We plot the dependence of gas phase metallicity 
(weighted by star formation rate), stellar mass, and gas 
mass, as function of circular velocity in Fig.[7]for various red- 
shifts. As expected, systems residing in shallower potential 
wells exhibit lower metallicities at all z, which is consistent 
with these galaxies being less able to retain metals in the 
presence of the kinetic energy injected by SNe in our model 
(see also Section roii . Interestingly the gas mass at a given 
circular velocity, does not evolve with 2. As a consequence, 
galaxies with given circular velocity have significantly higher 
specific star formation rates, and significantly lower Z, with 
increasing redshift. 

We have seen that the value of Z at a given stellar mass 
does not depend strongly on redshift. On the other hand, 
at a given Wire, galaxies tend to have significantly higher 
metallicities at lower 2. The evolution of the Z-V c i rc relation 
ranges from ~ 0.3 dex at low velocities to 0.6 dex or more 
at the high-velocity end. At first order, the evolution of the 
Z-V c irc relation can be understood as a consequence of the 
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Figure 7. SFR-weighted gas-phase metallicity (upper panel), 
stellar mass (middle panel) and gas mass (lower panel) as a func¬ 
tion V c i ic at various redshifts indicated in the legend. Curves de¬ 
note median relations, with error bars depicting the 25th and 
75th percentiles. At a given Wire, galaxies tend to have higher 
metallicities and higher stellar masses, at lower z, but similar gas 
masses. 


absence of significant evolution of the Z-M * relation but the 
significa nt evolutio n of t he M, -V c j rc relation (middle panel; 
see also ISwinbank et al.ll20l3 ). Galaxies with similar Wire 
(i.e. similarly deep potential wells) have a more dominant 
stellar component towards z ~ 0 and as a consequence of 
the MZR, higher metallicities. Because of the more efficient 
gas outflows of metals at low masses, the star formation 
in these systems is from less-chemical enriched material on 
average. We also showed that at a high Wire, the evolution of 
M* is stronger than in shallower potential wells. This issue 
can be probably related to the increase of the merger events 
for larger galaxies which do not only add the mass of the 
merging satellites but also contribute to the enhancement of 
the SF rate. 

In the lower panel of Fig. [7] we show the gas-phase 
mass in simulated galaxies as a function of Wire; interest¬ 
ingly there is almost negligible evolution. The increase of 


M* at a fixed potential well, with no significant variations 
of the gas mass, implies that the gas fractions of simulated 
galaxies decrease towards lower redshifts, consistent with the 
observed behaviour. This means that if the potential well of 
an individual galaxy does not change significantly over a 
given period of time and if this galaxy does not experience 
mergers over that period, its metallicity will increase as gas 
is converted into stars which enrich the ISM. In addition, as 
star formation proceeds consuming the gas-phase, gas infall 
will be required to maintain an almost constant gas compo¬ 
nent. 

4.3 Abundance imprints of galactic outflows 

As we have seen, our simulations are consistent with a sce¬ 
nario in which galactic outflows play an important role in 
setting the level of depletion of metals. We further investi¬ 
gate this process here by contrasting the abundances of el¬ 
ements produced by different channels, namely SNII, SNIa 
and AGB stars. Outflows in gimic are driven by SNe ex¬ 
plosions, which are assumed to occur ~ 30 Myr after star 
formation. Every star particle that is 30 Myr old imparts 
kinetic energy in the form of a randomly directed kick with 
a launch velocity of 600 km s _1 to, on average, r/ = 4 of its 
neighbouring gas particles. Type II, type I and AGB nucleo- 
synthetic produce are donated to neighbouring gas particles 
during each time-step. Since type II produce is distributed 
to neighbouring particles just before they are kicked, these 
elements are more likely to escape a galaxy in a galactic 
wind than elements produced by type I SNe and AGB stars, 
whose metals are released from much older stars and after 
the initial SNII kick. 

IWiersma et al- fl2010 ) discusses the enrichment of cos¬ 
mic metals in the OWLS suite, which used the same code, 
has comparable resolution to GIMIC for the high resolution 
runs, but were performed in periodic boxes. These authors 
find that metals residing in lower density gas were typically 
ejected earlie r an d by winds driven from lower mass galax¬ 
ies. IWiersma. Schave fc Theunsl (2011) investigate to what 
extent metal enrichment depends on the relatively poorly 
constrained sub-grid modelling of feedback. Their Fig. 12 
shows that in most models stars are the dominant reposi¬ 
tory of metals at redshift z = 0 (~ 60 — 70 per cent for the 
OWLS reference and MILL models for which the feedback 
implementation is most similar to that of gimic). However in 
the AGN model that includes feedback from accretion black 
holes, this decreases to 20 per cent. The second main repos¬ 
itory of metals in these models is the warm-hot intergalac- 
tic medium (WHIM). These results demonstrate that the 
metallicity of stars or star forming gas in GIMIC is not just 
set by the star formation efficiency or yields: metal mass loss 
resulting from feedback indeed plays an important role. For 
halos of mass below ~ 10 12 Mq feedback is regulated mostly 
by SNe driven winds, whereas in more massive haloes AGN 
feedback may be more important. Note that GIMIC does not 
include AGN. 

We plot the SF gas abundance X/Fe (X=N, Si, O, Mg, 
C and Fe from type I SNe only) as a function of M* and 
as function of Wire at z = 0 in Fig. [ 8 ] (left and right panel, 
respectively); abundances are off-set vertically so that the 
lines cross at M* = 1 O 1 O M 0 (Wire = 150 km s _1 for the 
right panel). The ratio a/Fe of a elements (a=0, Si, Mg 
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Figure 8. Median abundance of X/Fe for various elements X (X=N, Si, O, Mg, C and Fe from type la SNe only, as indicated in the 
legend), in star forming gas at redshift 2 = 0 as function of stellar mass M* and circular velocity Wire (left and right panels, respectively); 
error bars denote the 25th and 75th percentiles. For the sake of clarity, not all error bars are shown and we applied a small offset in the 
horizontal direction between some error bars to improve visibility. Abundances are off-set vertically by a constant W chosen such that 
they are all equal at M* = 10 10 Mq (Wire = 150 km -1 in the right panel), W = -0.685 (C), 0.127 (Mg), -1.08 (0), -0.14 (Si), 0.03 (N) 
and 0.36 (Fe from SNIa only). The ratio N/Fe increases dramatically with increasing M* or Wire, whereas the contribution of Fe from 
type I SNe decreases; these two elements bracket the behaviour of the others. 


which are produced in type II SNe) compared to (total) Fe 
increases with increasing mass up to M ~ 10 10 Mq, then 
flattens for O and Mg but keeps rising for Si; a similar trend 
is seen in a/Fe as function of Wire- This ratio depends on 
the star formation time scale, the delay time distribution of 
type I SNe, and whether or not type II produce is more likely 
to escape the galaxy in a wind than type I produce. However 
it of course also depends on nucleo-synthetic yields, which 
themselves depend on metallicity and hence stellar mass - 
this is the reason these elements do not track each other 
in detail in GIMIC . The X/Fe-M* relations for the stellar 
component are similar to those of star forming gas - modulo 
minor changes in the zero point. 

We mentioned in Section rm that simulated galaxies 
have around the same mass-weighted stellar age (~ lOGyr). 
Given the similar mean ages of our galaxies, the slope of the 
X/Fe-M* relation will be significantly influenced by the ef¬ 
fects of galactic winds at different masses. From Fig. [8] left 
panel, it is clear that there is a correlation between M* and 
X/Fe for C, Mg, O, Si and N and an anti-correlation between 
M* and the fraction of Fe from SNIa (/Fe.SNia). These results 
are consistent with what we would expect if low-mass sys¬ 
tems were more affected by galactic outflows. The shallower 
potential wells where low-mass systems reside are less effi¬ 
cient in retaining metals against galactic winds generated by 
SNII. On the other hand, due to the longer time scales for the 
release of SNIa yields, Fe from SNIa is not directly affected 
by these winds, generating lower X/Fe. It is worth noting 
that C and N are mainly generated by AG B stars b u t the 
contr ibution from SNII is also significant dwiersma et al.l 
120091) so that they are also affected by winds. In the right 
panel of Fig. [5] we can see the dependence of X/Fe on the 
potential well of simulated galaxies. Median X/Fe (/Fe.SNia) 
seem to be constant at log(W c i r c) < 2.1 and increase (de¬ 


crease) for larger log(W c i r c), indicating that galactic winds 
dominates at low Wire values, as expected. 

We note that the slope of the X/Fe-M* relation is shal¬ 
lower for C/Fe, probably caused by the significant contri¬ 
bution of AGB stars to C enrichment. A steeper slope is 
obtained for Si/Fe and N/Fe. I n the case of Si, this may be 
relate d to the assumed yields of lPortinari, Chiosi fe Bressanl 
d 1998t) for which Si production decreases with the metallicity 
of the progenitor star. 

The interpretation of the trends is more convoluted for 
the case of N which has characteristics of both a primary and 
a secondary elements. At low Z, N production results mainly 
from burning the C and O generated during the star’s evo¬ 
lution. Consequently its yield tracks the yields of C and O, 
and N behaves as a primary element However at higher Z, 
N production results increasingly from burning the C and O 
that was already in the star from birth and N behaves more 
as a secondary element - the N yield no longer tracks the 
C/O yield but depends on Z. Therefore, for galaxies with 
low metallicity, N is expected to behave as a primary ele¬ 
ment showing no important variations in N/O while for more 
metal-enriched galaxies, N/O should increase with O/H as 
N behaves as a secondary element. Because of the tightness 
of the MZR (Fig(4|, the abundance of O/H is almost a proxy 
for M*. This reasoning explains why the dependence of N 
abundance on M* steepens increasingly with increasing M*, 
as seen in Fig. [8] The behaviour of N/O with stellar mass 
may be different if Nitrogen were produced a bundantly in 
low-Z fast-rotating stars, as in the models of lPipino et all 
J2009I) . 

The study of the stellar X/Fe as a function of M* for 
elliptical galaxies is discussed extensively in the literature 
(see e.g., IConrov, Graves fc van Dokkumll2014l) . Although 
a detailed analysis is beyond the scope of the present paper, 
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Figure 9. Median mass fraction of Fe from SNIa in star forming gas as a function of M* {left panel) and Vf.irc {right panel) for various 
redshifts (different curves); error bars denote 25th and 75th percentiles. The contribution to total Fe from type I SNe increases with time 
as expected from the delayed production of type I Fe in the model. The dependence with M * (or V c ; rc ) due to more significant loss of 
SNII produce from shallower potentials, also becomes more pronounce at later times. 


we did perform a preliminary analysis of the stellar X/Fe- 
M* relation using a subsample of galaxies with large velocity 
dispersion, typical of elliptical systems. We find larger aver¬ 
age va lues of [M g/Fe], [S i /Fel a nd [N/Fe] than reported by 
IConrov. Graves fe van Dokkuml d20l4) by ~ 0.1 dex, ~ 0.3 
dex and ~ 0.15 dex, respectively. O bservatio nal u ncertain¬ 
ties a re of order 0.05 dex (IConrov, Graves fe van Dokkuml 
l20l4l and differences with other observational works are of 
order ~ 0.1 dex. Uncertainties in yields, and in the rates of 
SNIa, may be enough to explain the remaining differences. 

Finally we compare the abundance evolution of yields 
from SNII and SNIa as a function of Af* and I4irc (Fig. [9}, 
by plotting the fraction of Fe produced by SNIa. We see a 
significant trend of increasing SNIa Fe/total Fe with decreas¬ 
ing z, and decreasing M* (or Wire). The evolution with z is 
stronger for lower mass galaxies. These trends are consistent 
with our findings that type II products escape more easily 
from lower mass galaxies, and with the delay of SNIa enrich¬ 
ment. We reiterate that we may be overestimating the SNIa 
contributions especially for lower mass galaxies, because the 
stars in our low mass galaxies tend to be too old on average. 


5 SATELLITE GALAXIES 

In Section [5] we defined a central galaxy as the most mas¬ 
sive subfind structure in a given friends-of-friends halo - all 
other galaxies in that halo are satellites. Here we investigate 
abundance trends in satellites using the same criteria as for 
centrals regarding number of particles and mass range (see 
Section [2(2} . Consequently our sample will only include the 
most massive satellites; typically only 20 — 30% of the total 
depending on redshift. 

The MZR of satellites and central galaxies at z = 0 are 
compared in Fig. [TO] At given A/*, satellites are on average 
more metal rich than centrals by ~ 0.1 — 0.2 dex, and their 
metallicity scatter is higher than for centrals (~ 0.5 dex for 
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Figure 10. SFR-weighted gas-phase metallicity as a function of 
stellar mass at z = 0 for central (black) and satellite (red) galaxies; 
curves depict the median relations with error bars indicating the 
25th and 75th percentiles. At a given M*, satellites have higher 
median metallicities by ~ 0.1 — 0.2 dex, and their scatter in Z is 
also significantly higher. Both contribute to making the scatter in 
Z at given M* much higher for a sample of galaxies that includes 
both central and satellites compared to a sample consisting of 
just centrals. 


satellites compared to ~ 0.1 —0.2 dex for centrals), mainly at 
intermediate masses M* ~ 10 10 Mq. As we noted in Section 
EH the inclusion of satellites in the analysis of the MZR 
based on stellar abundances tends to increase the scatter. 
Our findings indicate that this is also true for gas abun- 
dan ces; t heir inclu sion also resu lts in an of fset to higher Z. 

IPasauali. Gallazzi fe van den Boschl (120121 ) compare 
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Figure 11. Median SFR-weighted gas-phase metallicity Z as a 
function of stellar mass for satellite galaxies at different redshifts 
2 ; error bars depict the 25th and 75th percentiles. Metallicity 
increases with M*, but unlike the case of centrals plotted in Fig. [5] 
Z also evolves with 2 , with satellites ~ 0.2 — 0.4 dex more metal 
rich at 2 = 0 compared to 2 = 3 at given stellar mass. The scatter 
in Z at given M* is also significantly larger than for centrals. 


the MZR of satellite and central galaxies inferred from SDSS 
data, finding that satellites tend to have higher metallicities 
than central systems at a given M* by ~ 0.06 dex, with 
the offset increasing with decreasing M*. We find a slightly 
larger offset, and no obvious trend with M* except at the 
very massive end. Previously we explained why varying the 
size of the aperture within which Z is measured affects the 
zero point of the MZR. because of abundance gradients - we 
find the same to be true for satellites. It is possible that this 
affects the interpretation of the data from SDSS due to the 
size of the fibres. 

What is the origin of the differences in Z between satel¬ 
lites and centrals? At a given M», satellites reside in shal¬ 
lower potential wells than central galaxies as expected: their 
Vcirc are lower by ~ 0.05 dex than those of central galax¬ 
ies of similar masses. Moreover, the / g -M* relation of satel¬ 
lite galaxies is offset by -0.2 when compared to centrals, 
mainly due to the loss of non star-forming gas. This lower 
gas fraction of satellites is consistent with these systems be¬ 
ing subject to strangulation and ram-pressure stripping. We 
also find that the scatter in / g -M* (and / g — Kirc) is larger 
for satellites than for centrals. These findings are consistent 
with a scenario in which satellite galaxies loose the metal 
poor gas in their outskirts due to ram-pressure stripping. 
This gas is then no longer available to dilute the ISM of the 
satellite, which leads to them having a higher Z. The larger 
scatter in / g is also consistent with the larger scatter in Z 
at given M* in satellites compared to centrals. Such a sce¬ 
nario can potentially explain both the offset and increase in 
scatter in Z, but we postpone a more detailed investigation 
to future work. 

The evolution of the MZR for satellites is shown in 
Fig. 1111 which can be compared to the evolution for cen¬ 
trals in Fig. [5] We recall that the ratio O/H - a good proxy 


for Z - does not show any appreciable evolution for centrals. 
In contrast Z does evolve for satellites, increasing from 2 = 3 
to 2 = 0 by ~ 0.2 — 0.4 dex, depending on M*. Interestingly 
at 2 ~ 3, satellites and centrals have MZRs with similar am¬ 
plitude and slope. We have verified that simulated satellites 
inhabit halos typically 10 times more massive than centrals 
at a given M*, which contributes to the higher metallicity 
of satellites. 


6 DISCUSSION 

Galaxies in the GIMIC simulation follow a mass-metallicity 
relation consistent with that observed, although the evolu¬ 
tion of the MZR may be less than that observed. We argued 
this maybe at least partly a consequence of bias in the ob¬ 
servations. As shown in Sections roi and R31 the slope of 
the MZR is to first order set by the efficiency with which 
galactic winds can escape the potential well. I 11 particular 
type II SN products such as Oxygen escape more readily 
from shallower potential wells. Since Oxygen is often used 
as a proxy for metallicity Z , ‘metals’ also escape more eas¬ 
ily from shallower potential wells. We note that this does 
not immediately imply that it is the escape velocity that 
determines whether metal-enriched winds can escape from 
a galaxy: higher cooli ng loses in a denser IS M m ay pla y an 
important role as well JCreasev, Theuns fe Bowerll2015l') . We 
also noted in Section nm that the simulated MZR does not 
exhibit the observed flattening at the high mass end, which 
we counted on the absence of AGN feedback in gimic. We 
conclude that feedback from SNe modulates the slope of the 
MZR at low and intermediate masses while the flattening at 
higher masses is increasingly affected by AGN. 

While galactic winds appear to be the main process for 
setting the slope of the MZR, its zero point is also affected 
by the SF time scales. The zero point of the X/H-M* rela¬ 
tion for the a elements produced by short-lived stars does 
not exhibit significant evolution since 2 ~ 3. We may have 
underestimated the level of evolution in a elements because 
our galaxies tend to be too old. However, the ratio X/H-A4* 
at given M* does evolve for elements produced by SNIs and 
AGB stars, and the level of enrichment is connected to the 
SF time scale. 

Applying observed aperture corrections and including 
satellites in our analysis increase the level of evolution by 
up to ~ 0.4 dex, however it is still not as strong as many 
observers infer. We believe this is at least partly due to the 
fact that stellar populations in the simulation are typically 
older than in the real universe, plausibly because feedback 
in the si mulation is not sufficiently effective at earl y times 
(see also lWeinmann et all201^ : lDe Rossi et aPl2013l l. In ad¬ 
dition we find that typical specific star formation rates in 
GIMIC are too low at redshift 2 = 0, in particu l ar for low- 
mass galaxies, as also found bv lMcCarthv et~ahl (12012) ') who 
blamed this partly on lack of numerical resolution. We note, 
however, that the overall MZR is not strongly affected by 
resolution, as demo nstrated in th e Appendix. In the recent 
EAGLE simulations (ISchave et alj 120151 ). feedback becomes 
more efficient at higher density which suppresses star for¬ 
mation at early times and increases the specific star forma¬ 
tion at later times. This may help to increase the level of 
evolution in the simulations. 
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With respect to the scatter in the MZR, our findings 
indicate that at a given M„, metallicities show a secondary 
dependence on / g , where galaxies with lower metallicities 
tend to have higher / g and higher SFRs, consistent with the 
observed trends. Therefore, the SF process in these systems 
is driven mainly by accretion of metal poor gas. 

Interestingly, we find that, at a given M*, systems with 
lower gas-phase metallicities tend to exhibit higher Kirc- 
Therefore, at a given M*, systems with deeper potential 
wells have larger fractions of gas which is also more metal- 
poor. Such a correlation suggests that the depth of the po¬ 
tential well at a fixed M* is regulating the infall of metal- 
poor gas and hence, the scatter in the MZR. In addition, 
at a fixed M*, the radius which encloses 50% of the stellar 
component (R 50 ) tends to be larger for systems with lower 
metallicities, evidencing again the more significant influence 
of gas infall in these galaxies. 

In summary: the zero point of the MZR in the GIMIC 
simulations is at first order set by the star formation time 
scale, and the efficiency of galactic winds, respectively. The 
scatter around the relation is affected by the gas fraction, 
itself correlated with the depth of the potential of a galaxy 
of given mass. The MZR of satellites is off-set to higher Z 
compared to centrals of the same mass, and its scatter and 
level of evolution are significantly higher. We verified that 
ram-pressure stripping and the metal pollution of the host 
halos in which satellites live, seem to play a key role on 
modulating their metal enrichment. 


7 CONCLUSIONS 

We have analysed the GIMIC hydrodynamical simulations of 
structure formation to study the abundances of different ele¬ 
ments in star forming gas in galaxies as a function of galaxy 
stellar mass and redshift. We focused mainly on the anal¬ 
ysis of central substructures of dark matter halos (‘central 
galaxies’) but we also looked at satellite systems. Our main 
conclusions can be summarised as follows: 

• GIMIC galaxies exhibit a well-defined correlation be¬ 
tween O/H and M» at 2 = 0 in the sense that more mas¬ 
sive galaxies are more metal-enriched, in agreement with 
observations (Fig. [3. However, at the high-mass end (M* > 
10 10 - 5 Mq) simulations predict higher abundances than ob¬ 
servations. We suggest that the inclusion of AGN feedback 
in the model would help to overcome this problem by pre¬ 
venting overcooling in massive galaxies. 

• Consistent with observations, the scatter of the O/H- 
M, relation is driven by a secondary dependence of metal- 
licity on the gas fraction (/ g , Fig. [2]). At given M*, galax¬ 
ies with lower metallicities Z tend to have higher / g , higher 
star formation rates, inhabit deeper potential wells and have 
larger half-mass radii. Such correlations are expected if star 
formation in galaxies is self-regulating (ISchave et al.l [20101 1 
with feedback regulating inflow of more metal poor gas. 

• The simulations are able to reproduce the correla¬ 
tion between O/H, M* and SFR of the fundamental 
mass-metallicity relat ion (FMR) of observed galaxies of 
iMannucci et al l d2Qloi f ■ Our results are also consistent with 
the existence of a correlation between O/H, M* and / gas . 
Following observational parametrizations of the FMR, we 
analysed the parameters fj, a and rig. We found that /3 ~ 0.48 


leads to the tightest relation with O/H, reducing the scatter 
by ~ 0.03 dex with respect to the O/H-M* relationship. Our 
findings show also evidence for a correlation b etween / g , 
and S FR, consistent with the observations bv ISantini et al.l 
J2014ll . At a given M*, simulated systems with higher / g 
tend to have higher SFRs. The simulated fg-M * relation for 
the cold gas is in remarkably good agreement with observa¬ 
tions. 

• The Z*-M* relation in GIMIC agrees well with obser¬ 
vations in both slope and zero point (Fig. []}; however the 
scatter in the simulations is lower. Including satellites in the 
analysis, and taking into account aperture effects increases 
the scatter in the simulation considerably. 

• There is negligible evolution of the O/H-M* relation 
between 2 = 3 and 2 = 0 (Fig. [SJ), in apparent disagreement 
with observations. Taking into account satellite galaxies - 
that do evolve - and aperture effects leads combined to a 
small evolution of ~ 0.4 dex from 2 = 3 to 2 = 0, still 
considerably less than the claimed 1 dex evolution in the 
data. In GIMIC, galaxy stellar ages tend to be uniformly old, 
~ 10 Gyr which contributes to the disagreement - modi¬ 
fying the ste llar feedba ck efficie ncy as in the recent EAGLE 
simulations dSchave et al.ll2015l'l might resolve some of the 
discrepancy. However we suggest that observational bias in 
the data may be the main culprit. 

• The evolution of the abundance ratios (X/H) as a func¬ 
tion of galaxy stellar mass, for various other elements X 
tracked in the simulation, depends strongly on the stel¬ 
lar evolutionary channel that produces X (Fig. [5J. a- 
elements produced by massive stars show very little evolu¬ 
tion, whereas elements with significant SNI or AGB contri¬ 
butions show stronger evolution. Therefore, the total metal- 
licity (Z) of the gas evolves in the sense that at a given 
mass, systems were less enriched in the past. The different 
behaviour of these elements can be understood in terms of 
(*) the ejection of SNII produce in galactic winds, and ( ii ) 
the stellar evolutionary time scale for SNI and AGB stars. 
Galactic winds are driven by the same type II SNe that also 
release a elements, and such metal enriched winds escape 
easily from small galaxies living in a shallow potential well. 
If the stars in these galaxies are old enough to produce SNI 
and AGB stars, then elements produced predominantly by 
SNI and AGB stars dominate the abundance pattern. More 
massive galaxies from which winds cannot escape easily are 
then relatively more abundant in a elements. We focused 
mainly on the comparison of the slope of the different abun¬ 
dance ratios as the theoretical yields are uncertain by a fac¬ 
tor of a few. A comparison between results presented here 
and those derived from other yields models available in the 
literature would be very interesting. 

• Satellites are more metal rich than centrals at given M* 
in agreement with observations, and the scatter around the 
mean relation is significantly larger for satellites ('Fig. 1101) . 
The inclusion of satellites in the analysis of the MZR could 
increase the scatter of the relation, leads to an offset, and 
introduces evolution. 

Satellites have similar amounts of star forming gas, and 
similar star formation rates as central galaxies at a given 
M*. However, they have less gas in total, probably be- 
cause of strangu l ation and ram-pressure stripping (see also 
iBahe et aLll2013l : lBahe fe McCarthvll2015l . for a more exten- 
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sive study of the satellite galaxy population in GiMic). Their 
depletion in gas contributes to their higher metallicities. 

The mass-metallicity relation in GIMIC is mostly deter¬ 
mined by the star formation time scales and the efficiency of 
SN feedback as a function of galaxy stellar mass. SN feed¬ 
back is efficient in shaping the relation at low and inter¬ 
mediate masses, A/* < 10 10 ' 5 M@, but at higher masses the 
absence of AGN feedback in the simulation leads to unreal¬ 
istically massive galaxies. The evolution of the zero point of 
the MZR seems to be directly related to the SF time scales. 
Since GIMIC galaxies have a dominant older stellar compo¬ 
nent, the MZR based on a elements alone does not evolve 
significantly. We suggest that an improvemen t of t he feed¬ 
b ack prescription as implemented in EAGLE dSchave et al.l 
l2015ll might affect this. At second order, the simulated MZR 
shows dependencies on the infall of metal-poor gas. At a 
given M*, systems with deeper potential wells seem to acrete 
more metal-poor gas increasing their SFRs. These findings 
are consistent with previous proposed scenarios in the liter¬ 
ature. 

Finally, although there are still important discrepancies 
between observational studies with respect to the shape, zero 
point and evolution of the MZR and its 3D extensions, the 
advances of the observational techniques are shedding more 
light on the abundance evolution of galaxies over cosmic 
time. The continued comparison between these observations 
with theoretical models of galaxy formation not only offers a 
fruitful means for testing/constraining the models, but also 
helps to guide those studies and provide plausible formation 
scenarios. A significant combined effort will be required from 
the theoretical and observational sides in order to converge 
to a complete understanding of the abundance enrichment 
history of galaxies. 
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APPENDIX: RESOLUTION STUDY 

The analysis presented in this paper has been done by using 
the high resolution version of GIMIC simulations (see Section 
H. In order to assess the numerical convergence of our re¬ 
sults, we have also performed a similar analysis by using the 
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Figure 12. SFR-weighted gas-phase metallicity as a function of stellar mass at 2 : ^ 0 for central galaxies in the high resolution run 
(black) and intermediate resolution run (red). The different curves with error bars depict the median relation with the 15th and 85th 
percentiles. The level of numerical convergence is good. 
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Figure 13. Stellar metallicity as a function of stellar mass at z ^ 0 for central galaxies in the high resolution run (black) and intermediate 
resolution run (red). The different curves with error bars depict the median relation with the 15th and 85th percentiles. Convergence is 
good, with a small tendency for higher values of Z at higher resolution. 


intermediate resolution runs, which has 8 times coarser mass 
resolution. Our results are robust to a factor of 8 change in 
mass resolution 

As an example, we compare the MZR based on SF- 
weighted gas and stellar metallicities at two resolutions, in 
Figs. I12l and ll.il respectively. We see that the main features 
of the relations seem to be robust against resolution. There 
is only a very small trend in the case of the intermediate 
resolution run to have lower abundances at the high-mass 
end towards higher redsliifts. 
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